!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! the subroutine is used to calcaulate direct wall 
! papramenters by MOC.
! ref. to chapter 16 of 'Gas Dynamics' by Maurice J. Zucrow & Joe D. Hoffmann
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! coded by : B. G.
! created  : 2015-06-17
! revised  : 2015-07-06 add cubic equation for wall function
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! include  :
!   subroutine : DirectW
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! nomenclature:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! warning:
!   1. angle in rad, not degree
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine DirectW()
   use VariableDef
   implicit none
!     (1)
!            4()
!     2
!     given variables
!   real*8::xa,ya,thetaa,thetae,xe
!   real*8::x2,y2,Vu2,Vv2
!     temp variables
   real*8::thetap,ap,alp,Vup,Vvp,Vp,yp,Mp
   real*8::lp,Qp,Rp,Sp,Tp
   real*8::a4,x4f,y4f,Vu4f,Vv4f
   integer::icor,i
!,wallflag
!     calculate variables
!   real*8::x4,y4,Vu4,Vv4,V4,p4,T4,rou4,M4
   
   x4f=0.0
   y4f=0.0
   Vu4f=0.0
   Vv4f=0.0
   icor=6
   ya=abs(ya)
   thetaa=abs(thetaa)
   thetae=abs(thetae)
   y2=abs(y2)
   Vv2=abs(Vv2)

   Vup=Vu2
   Vvp=Vv2
   yp=y2
      
   do i=1,icor
      thetap=atan(Vvp/Vup)
      Vp=sqrt(Vup**2+Vvp**2)
      call Thermo(Vp,p4,rou4,T4,ap,Mp)
      alp=asin(1/Mp)
      
      lp=tan(thetap+alp)
      Qp=Vup**2-ap**2
      Rp=2*Vup*Vvp-lp*Qp
      Sp=delta*Vvp*ap**2/yp

      call Wall24(lp)   ! calculate 4 location by second order wall funtion
!      call Wall34(lp)   ! calculate 4 location by third order wall funtion
      
      Tp=Sp*(x4-x2)+Qp*Vu2+Rp*Vv2

      Vu4=Tp/(Qp+bw*Rp+2*cw*Rp*x4)   ! velocity calculated by second order wall funtion
      Vv4=Tp*(bw+2*cw*x4)/(Qp+bw*Rp+2*cw*Rp*x4)
!      Vu4=Tp/(Qp+Rp*(bw+2*cw*x4+3*dw*x4**2))   ! velocity calculated by third order wall funtion
!      Vv4=Tp*(bw+2*cw*x4+3*dw*x4**2)/(Qp+Rp*(bw+2*cw*x4+3*dw*x4**2))
         
      if ((abs(x4-x4f) <= el) .or. (abs(y4-y4f) <= el)) then
         if ((abs(Vu4-Vu4f) <= ev) .and. (abs(Vv4-Vv4f) <= ev)) then
            exit
         end if
      end if

      Vup=(Vu2+Vu4)/2
      Vvp=(Vv2+Vv4)/2
      yp=(y2+y4)/2
      
      x4f=x4
      y4f=y4
      Vu4f=Vu4
      Vv4f=Vv4

   end do
   V4=sqrt(Vu4**2+Vv4**2) !point 4 velocity
   call Thermo(V4,p4,rou4,T4,a4,M4)
      
   if (wallflag==-1) then
      ya=-ya
      thetaa=-thetaa
      thetae=-thetae
      y2=-y2
      Vv2=-Vv2
      y4=-y4
      Vv4=-Vv4
   end if
end subroutine DirectW
